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We develop a new numerical scheme which allows precise 
solution of coherent tunneling problems, i.e., problems with 
exponentially small transition amplitudes between quaside- 
generate states. We explain how this method works for the 
single-particle (tunneling in the double-well potential) and 
many-body systems (e.g., vacuum-to- vacuum transitions), 
and gives directly the instanton shape and tunneling ampli- 
tude. Most importantly, transition amplitudes may be calcu- 
lated to arbitrary accuracy (being limited solely by statistical 
errors) no matter how small are their absolute values. 
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Tunneling phenomena are among the most intriguing 
consequences of quantum theory. They are of fundamen- 
tal importance both for the high-energy and condensed 
matter physics, and the list of systems which behavior is 
governed by tunneling transitions ranges from quantum 
chromodynamics jjjfl to Josephson junctions and defects 
in crystals (see, e.g., Ref. ||). 

Precise analytic treatment of tunneling in complex sys- 
tems is very hard, if not impossible. The most crucial 
simplification is in reducing the original problem to the 
semiclassical study of the effective action for some collec- 
tive variable R, with the assumption that all the other 
degrees of freedom adjust adiabatically to the motion of 
R HH . In certain cases, one may also include dissipa- 
tive effects due to "slow" modes other than the selected 
collective variable which do not follow the dynamics of 
R adiabatically Typically, the parameters of the 

effective action can not be found analytically (although 
one may relate some of them to the linear response coeffi- 
cients) and have to be deduced from experiments. As far 
as we are aware, at present there are no tools to address 
the tunneling problem numerically. "Exact" diagonaliza- 
tion works only for relatively small systems, and, even in 
small systems, its accuracy is not sufficient to resolve very 
small energy splittings AE, say, when AE/E <C 10~ , 
due to round-off errors (unless specific no-round-off arith- 
metics is used). 

In this letter we develop a quantum Monte Carlo (MC) 
approach which allows precise calculations of tunneling 
amplitudes (and instanton shapes) no matter how small 
are their absolute values. The MC scheme contains no 
systematic errors, and its accuracy is limited only by sta- 
tistical noise. The key point of our approach is in simulat- 
ing imaginary-time dependence of transition amplitudes 
Aji(r) between selected reference states and 1 772 ) (see 
below) . In doing so we have to solve the problem of col- 



lecting reliable statistics in a case when Aji(r) varies, 
say, over hundreds of orders of magnitude (!) between 
different points in time. First, we describe in detail how 
to evaluate tunneling splitting and instanton shape in 
the double-well potential. We proceed then to the many- 
body problem of vacuum-to-vacuum transitions by con- 
sidering the case of ID quantum antiferromagnet with 
exchange anysotropy. Finally, we discuss the generality 
of the method suggested and demonstrate our numeric 
results for tunneling in the double-well potential (with 
comparison to the exact-diagonalization data where pos- 
sible) . 

Consider the standard problem of particle motion in 
external potential: 



H = mx 2 /2 + U(x) 



(1) 



= rji and x 



which has two minima at points x 
r]2, and large tunneling action S , 

p 2 [2mU(x)]^ 2 dx > 1. These minima are supposed 
to be near-degenerate, i.e., the lowest eigenstates of the 
Hamiltonian (111), H^ r 



£ a f „, form a doublet with 
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LJi < LU t 
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where are the classical vibration frequencies in the 
potential minima u>i = [U" (rj^mj 1 ^ 2 (we set fi = 1). The 
lattice analog of the Hamiltonian ([!]) reads 



H = -t d\d v + ^2 niUi , ni = djdi , 



(3) 



<W> 



where dj creates a particle on the site number I, and the 
first sum is over nearest-neighbor sites. 

The transition amplitude from the state \r/i) to the 
state |r72 ) , where ^1,2) = S(x — 771.2), in imaginary time 
r is given by 

Ajiir) = (m\e- HT \m) = J2( a \K)(V2\<x)e- E « T . (4) 

a 

We now make use of the inequality (||) to define the 
asymptotic regime E<x — E\ -C r _1 <C Wj, see, e.g., 
Refs. §, 

A ., {t ) _ e -Er J- (a\ m )( m \a)[l - (E a - E)r] . (5) 

a=l,2 

where E=(E 2 + E 1 )/2. 

It is convenient to split the double- well potential in two 
terms U(x) = U {1) {x) + U i - 2) {x), where U {1 ' 2 \x) is identi- 
cal to U (x) to the left /right of the barrier maximum point 



1 



and remains constant afterwards. Introducing system 
ground states in each minimum as H (l >^ G = E G ^ G , 
where H^"' — mx we may rewrite 



(i) 



it* 



(2) 
G ' 

(2) 



-G _ J ( 6 ) 

where (w 2 ,i; 2 ) = l/2±£/2E, E 2 = A 2 +£ 2 , with obvious 
identification of the energy splitting 2E = E 2 — E\ and 
bias energy i = (E ( G ] - E [ G ] )/2. Here A is the tunneling 
amplitude, which defines energy splitting E 2 — Ex = 2A 
in the degenerate case. Substituting Eq. (||) into Eq. (||) 
we finally obtain 



Au[t) 
A jM t ) 



-Et ryl 



(7) 



Here Zi 



W 



7/i) projects the reference states 1 77^ ) on 



the corresponding ground states in each minimum. All 
corrections to Eq. (^), e.g., the neglect of the overlap in- 



tegrals (*g ) |?72) and (* 



(2) 



771}, are small in parameter 



e 6 . Note also that Eq. (0) does not depend on the bias 
£ <C u>i, and is insensitive to the behavior of ^ G in the 
deep underbarrier region. This (rather standard) con- 
sideration relates transition amplitudes to the tunneling 
amplitude through the asymptotic analysis of A^ (r) in 
time. 

MC simulation of the transition amplitude is almost 
identical to the standard simulation of quantum statistics 
[the partition function at a given temperature T = 1/(3 
may be written as Z[r = /3) = Tt^j4(?7, rj, r)]. Fixed 
boundary conditions, as opposed to the trace over closed 
trajectories (configurations), are trivial to deal within 
any scheme. The real difference is in sampling different 
time-scales - thermodynamic calculations are typically 
done with (3 = const. Now we are forced to consider tra- 
jectories with different values of r and to treat imaginary 
time as "dynamic" (in MC sense) variable. The idea of 
utilizing time dependencies of trajectories in MC simu- 
lations was extensively discussed in connection with the 
Worm algorithm {jjJI and polaron Green function ||. 

We now turn to the problem of normalization. This 
problem might seem intractable in view of close analogy 
between Z(t) and A(t). Formally, in the limit r — > 0, 
the amplitude is trivial to find in most cases, e.g., for the 
Hamiltonians ([!]) and (||) it is given by the free particle 
propagation 

-%) 2 /2r 




continuous 
discrete , 



(8) 



and this knowledge may be used to normalize MC statis- 
tics for Aji(r) (in the discrete case \rji — rjj | =integer). 
However the absolute values of Aji (r) at short and long 
times will typically differ by orders and orders of magni- 
tude and simply none MC statistics will be available at 
short times. 



The solution to the puzzle lays in the possibility to 
use an arbitrary fictitious potential A^ c (t) in Metropolis- 
type updates pOl in time-domain 



P acc {B -> A) 
P acc (A -> B) 



-3 B Aic(r') W A (r) 
Aflc(r) W b (t') 



(9) 



where P acc (B — > A) / P aC c(A — > B) is the acceptance ratio 
for the update transforming initial trajectory B (having 
duration r) to the trajectory A (having duration r'), and 
g-Se j g fo e statistical weight of the trajectory B. The 
normalized distribution functions Ws(r'), according to 
which a new value of r is seeded, are also arbitrary; the 
best choice of W's follows from the conditions of (i) opti- 
mal acceptance ratio (as close to unity as possible), and 
(ii) simple analytic form allowing trivial solution of the 

equation JJ" W(t')(It' = r J 71 W(t')cIt' on the time in- 
terval (r a ,rb) , where < r < 1 is the random number 
. Each trajectory adds a contribution = 1/Ar c (t) to 
the time- histogram of A^(t). 

One may use fictitious potential to enhance statistics 
of trajectories with certain values of r "by hand", e.g., 
by making A^ c zero outside some time-window. To get 
a reliable and properly weighted statistics both at short 
and long times we need Ar c (t) ~ 1/Ay(r) to compen- 
sate completely for the severe variation of the transi- 
tion amplitude between different time-scales. This goal 
is achieved as follows. The initial stage of the calcula- 
tion, called thermolization, prepares the fictitious poten- 
tial using recursive self-adjusting scheme - starting from 
Afic(f) = 1 in a given time- window (r m i„, T max ) and zero 
otherwise, we collect statistics for Ay(r) to the tempo- 
rary time-histogram and after M > 10 6 -j- 10 7 updates we 
renew the fictitious potential as 



Aij (t )/ Aij (t) , 
Aflc(r) = <( A^ (t )/ A^ (n) , 
A^ (t )/ A^ (t 2 ) , 



T\ < T < T 2 
7~min ^ T <^ T\ 
T 2 < T < T max 



(10) 



where t\ and t 2 are the points (to the left and to 
the right of some reference point to) where temporary 
statistics becomes unreliable and has large fluctuations. 
It makes sense to select tq close to the maximum of 
Ajj(r) (this point may be tuned a posteriory), and points 
T\ and t 2 are formally defined as the first points in 
the histogram where smooth variation of A^(t) ends: 
Aj(ti,2 + Ar)/Aij(Ti t 2) < <5 (here Ar is the difference 
between the nearest points in time-histogram, and S is a 
small number, say 0.01). The thermolization stage con- 
tinues until Ti = T m i n and t 2 = T max , and fictitious po- 
tential stops changing (withing a factor of two). After 
that the actual calculation starts with a fixed A^ c (t), 
and a new histogram for Aij (r) is collected. 

The idea of using fictitious potential proportional to 
the inverse of the transition amplitude is clear - it allows 
to collect reliable statistics on different time-scales with 
comparable relative accuracy. With this tool at hand one 
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can easily normalize Ay (r) using known analytic results 
for short times, and deduce tunneling amplitude and Z- 
factors from the analysis of the long-time asymptotic, 
Eq. (@). We can not but note that fictitious potential in 
time-domain very much resembles the so-called "guiding 
wavefunction" in the Green- function MC methods [ Q |, 
with essential difference that here it is used to reach ex- 
ponentially rare configurations. 

Instanton shape calculation is a much easier task since 
it can be done by considering trajectories for the tran- 
sition amplitude Ai^j(r) with fixed but sufficiently long 
r 3> 1/w^; again, to ensure that Ay is dominated by 
just one instanton trajectory we need r <C 1/A. For any 
given MC trajectory x(r) one has first to define the in- 
stanton center position in time, and to recount all times 
from this center-point r c [instanton center statistics is al- 
most uniform in (0, r) (except near the ends of the time 
interval) due to the generic "zero-mode" present in in- 
stanton solutions [§,||]. This can be done by looking at 
the average time 



I B TdT 
I B dT 



(11) 



where the integral is taken over the barrier region be- 
tween the wells U{x(t)) > E G . The instanton shape 
is obtained then by collecting statistics of x(t — r c ) to 
the time histogram. In this simple example when the 
notion of collective coordinate is not necessary (or for- 
mally, R = a;), we do not need to define separately the 
estimator for R (see the opposite example below). 

The case of vacuum-to- vacuum transition in the many- 
body problem is formulated in precisely the same man- 
ner, and Eqs. (f||), (|-|o|) holds true once the identifica- 
tion of the states |i7i 2) is done and the short-time limit, 
Eq. (||), is calculated. Consider as a typical example a 
ID spin chain with 2L sites and antiferromagnetic (AF) 
couplings between the nearest-neighbor spins 



H 



E 

<ij> 



(12) 



which has near degenerate ground state with the lowest 
doublet separated from the rest of the system spectrum 
by finite gap for J' > 0. The natural choice of ^1,2) 
is then an ordered AF state with Sf |t?i,2) = (±^) l \Vi,2)- 
Note, that for a large system Z-factors \r)i) will be 
also exponentially small, and even diagonal amplitudes 
An have to be calculated with the use of A^ c [in the 
single-particle case one may obtain Z-factors by ignoring 
Ag c -trick] . The short-time behavior is given by 



Au(t) -» 1 
A^(t)^(Jt/2) J 



2 (ring) 

1 (open chain) 



(13) 



such knowledge is not available one has to study differ- 
ent possibilities). For example, if tunneling proceeds via 
two domain walls well-separated from each other (thin- 
wall approximation |i]j|||), then the collective variable 
R is the distance between the walls, and the underbar- 
rier region in Eq. (jll]) is related to the existence of two 
separated walls. These definitions are not too specific 
and work only approximately; this is however the generic 
difficulty of dealing with collective variables which are 
meaningful only in the macroscopic limit. Obviously, the 
knowledge of the instanton shape does not allow precise 
evaluation of A, and gives only rough estimation of In A. 

To test the proposed scheme and to compare results 
to the exact diagonalization (ED) data we have applied 
our algorithm to the lattice model (||) with U(x) = 
Uo[(x/r/) 2 — l] 2 . In what follows we measure all ener- 
gies in units of the hopping amplitude t and count them 
from the potential minimum. We set Uo = 1 and consider 
two inter well separations: 77 = 10 and rj — 40. For r\ = 10 
the ED data for the ground state energy, Z-factor, and 
tunneling amplitude are: Eq = (Ei + E 2 )/2 = 0.1923, 
Z 2 = 0.2465, A = (E 2 — E x )/2 = 3.6078 x 10~ 6 . 
Our MC data give E G = 0.192(2), Z 2 = 0.246(2), and 
A = 3.61(3) x 10 -6 . The case of large rj is more subtle 
since only Eq — 0.0495 and Z 2 = 0.1255 may be tested 
against ED - one may use textbook semiclassical anal- 
ysis of the corresponding continuous model to see that 
A (77 = 40) ~ 10~ 24 -j- 10~ 23 , that is far beyond the stan- 
dard computer round-off errors. 



Ln(Aii(T)) 




To decipher the instanton we need now some a pri- 
ori knowledge about the relevant collective variable (if 



FIG. 1. Time dependence of the diagonal amplitude An. 
Solid line Z 2 e~ E ° T is the prediction of the long-time behavior 
according to the ED data. 

In Fig. ([j]) and Fig. (J2J) we present our MC data for 
the diagonal and off-diagonal amplitudes for rj = 40, and 
fits to the expected long-time, Eq. (0), and short-time 
behavior, Eq. (H) [we also include the lowest-order cor- 
rection for the potential energy at r — > which tells that 
Aijij oc e _c/r where U = f_ U(x)dx\. The variation 
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of Ai£i in Fig. 2 is about two-hundred orders! From 
Fig. @ we deduced £ G = 0.0494(2) and Z 2 = 0.126(1) 
in agreement with ED. From Fig. (|2|) we then obtain 
A = 1.7(4) x 10 -23 . MC simulation of the above figures 
took about 5 days each on PII-266 processor. 
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the requirement of smooth variation of the amplitude 
Ardln[Ai^(r)]/dT < 1, after substituting Eq. (|), 
means At <C t/(2t]) ~ 10~ 3 for r = 0.1. To avoid 
large systematic errors due to time-discretization at short 
times one has to use Ar = 10~ 4 (!!). Apart from enor- 
mous memory usage (there will be about 10 7 time slices) 
that small Trotter parameter severely slows down the ef- 
ficiency of the code, in fact, we are not aware of any MC 
simulation with At ~ 10 . 

It is worth mentioning that similar technique makes it 
possible to study directly Z(j) over different time-scales 
- normalization of the partition function does not matter 
since none physical quantity depend on it. Thus one 
can obtain temperature dependences of the free energy, 
entropy etc. in a single MC run. 

This work was supported by the RFBR Grants 98-02- 
16262 and 97-02-16548 (Russia), IR-97-2124 (European 
Community). 



FIG. 2. Time dependence of the non-diagonal normal- 
ized amplitude A^j. The solid line is fitting to the curve 
AZ 2 e~ E ° T and the short-time behavior is tested against the 
law (r 2rl - 1 /(2r 1 )\)e' ijT (dashed curve). 
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FIG. 3. Instanton trajectory x(t). The solid curve 
x(t) = r\ tanh(_Ec r) is the semiclassical result for the con- 
tinuous model. 

In Fig. ([|) we present our results for the instanton tra- 
jectory x(t) (at r\ = 40) along with the known semiclassi- 
cal results for the continuous model (0) . The accuracy 
of the data is self-evident, although we argue that MC 
data rather represent the trajectory with finite energy 
E = Eg while analytic results correspond to E — 0. 

We note that the present technique is very hard 
to implement in the discrete-time schemes with fi- 
nite Trotter parameter At. On one hand, long-time 
asymptotic regime (see Fig. 2) requires to consider t 
as long as 640. On another hand at short times, 
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